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ABSTRACT 

The  present  paper  investigates  sources  of  uncertainties  in  two-dimensional  flow  compu¬ 
tations  and  presents  methods  for  estimating  them.  Two  sample  problems  are  used  for  illus¬ 
tration.  The  following  categories  are  explored  in  detail:  i.)  Uncertainty  due  to  truncation 
error  in  numerical  schemes;  ii.)  Uncertainty  due  to  discretization  error;  iii.)  Uncertainty  due 
to  outflow  boundary  conditions;  iv.)  Uncertainty  due  to  incomplete  iterative  convergence; 
V.)  Uncertainty  due  to  computational  grid  aspect  ratio.  The  error  estimates  are  based  on 
requirements  for  internal  consistencies  in  computed  results.  Therefore,  they  provide  bet¬ 
ter  judgement  of  the  numerical  solution  integrity  than  comparisons  to  experimental  data 
or  “benchmark”  solutions  whose  reliability  may  sometimes  be  questionable.  Ideally,  both 
approaches  should  be  employed.  New  results  are  presented  on  the  optimum  grid-cell  aspect 
ratio  for  computational  accuracy  and  efficiency. 

one  QUALIIT  INCFECTED  8 


Accesioii  For  | 

NTIS  CRA&I  M 

OTIC  TAB  B 

Unannouiiced  □ 

Justification 

B) 

Di 

f 

sti  ibution  / 

Availability  Codes 

Dist 

tL 

Avail  £ 
Spe 

ind/or 

cial 

'Research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract 
Nos.  NASl-18605  and  NASl-19480  while  the  first  author  was  in  residence  at  the  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681- 
0001. 


I 


1  INTRODUCTION 


The  rapid  development  of  computers  over  the  past  three  decades  has  encouraged  the  development 
of  computational  fluid  dynamics  to  such  an  extent  that  it  has  become  a  viable  analytical  tool  in  the 
solution  or  design  process  in  several  engineering  and  environmental  applications.  As  the  investi¬ 
gated  flow  situations  have  become  more  complicated,  the  need  for  techniques  for  evaluating 
sources  and  magnitudes  of  uncertainties  in  computed  results  have  grown.  Unfortunately,  too  little 
attention  has  been  paid  to  this  subject  in  the  literature.  Most  work  is  done  by  numerical  analysts  in 
highly  idealized  situations.  The  dilemma  facing  the  computational  fluid  dynamists  is  that  the  log¬ 
ical  test  of  the  accuracy  of  a  numerical  method  is  to  compare  computed  results  to  exact  solutions. 
But  these  can  rarely  be  found  for  practically  interesting  problems.  So  one  resorts  to  comparing 
computed  results  to  other  numerical  solutions,  sometimes  called  "benchmark"  solutions,  or  exper¬ 
imental  data.  Such  a  comparison,  may  be  diagnostic,  but  would  usually  be  inconclusive  since  the 
latter  may  contain  unknown  errors.  Therefore,  there  is  a  need  for  techniques  which  allow  quanti¬ 
tative  estimation  of  various  uncertainties  in  the  numerical  solution  in  a  systematic  manner. 

Ferziger  (1989)  proposed  some  methods  suitable  for  the  estimation  and  reduction  of  numerical 
errors  resulting  from  inadequate  grid  resolution  or  incomplete  convergence  of  the  iterative 
scheme.  The  former  is  based  on  the  Richardson  extrapolation  method  originally  proposed  by 
Richardson  (1911)  and  Richardson  and  Gaunt  (1927).  This  method  has  been  used  in  a  wide  range 
of  applications  to  improve  numerical  solutions  or  to  estimate  errors  in  numerical  solutions. 
Churchill  et  al.  (1981)  and  de  Vahl  Davis  (1983)  applied  the  method  to  estimate  zero-grid-size 
solution  in  natural  convection  problems.  Applications  to  aerodynamic  flows  are  reported  by  Dang 
et  al.  (1989)  and  Zing  (1991),  amongst  others.  The  common  result  is  that  the  Richardson  extrapo¬ 
lation  method  is  reliable  only  when  the  numerical  solutions  on  the  different  grids  used  in  the  pro¬ 
cedure  are  smooth  and  display  similar  characteristics,  which  presupposes  that  the  grids  should  be 
sufficiently  fine  to  resolve  all  flow  features.  The  method  for  removing  the  uncertainty  of  incom¬ 
plete  convergence  of  the  iterative  scheme  is  to  base  the  convergence  criterion  not  on  the  change  in 
computed  results  between  iterations,  as  is  common  in  practice,  but  on  an  estimate  of  the  solution 
error  constructed  from  both  the  change  and  its  rate  of  change.  This  method  requires  little  addi¬ 
tional  effort  but  has  not  yet  found  wide  use. 

In  the  present  paper,  we  investigate  a  wider  range  of  sources  of  uncertainty  in  numerical  computa¬ 
tions  of  separated  flows.  Possible  errors  resulting  from  each  source  are  estimated  and  methods  for 
eliminating  or  minimizing  them  are  explored.  Two  model  problems  were  used  in  the  investiga¬ 
tion.  The  first  is  the  steady,  two-dimensional  laminar  flow  over  a  backward  facing  step,  and  the 
second  is  a  steady,  two-dimensional  stratified  laminar  flow  over  a  backward  facing  step.  The  Rey¬ 
nolds  number  in  both  flow  problems  was  at  the  high  end  of  the  laminar  flow  range  (equal  to  400 
based  on  the  mean  flow  velocity  and  the  channel  height  upstream  of  the  step).  The  flow  configura¬ 
tions  and  the  boundary  conditions  are  illustrated  in  Figures  1  and  2.  “Benchmark”  solutions  for 
these  problems  have  been  developed  by  Gartling  (1990)  for  problem  1  and  Leone  (1990)  for  prob¬ 
lem  2.  Streamlines  of  these  are  shown  in  Figure  3.  Our  results  are  also  compared  to  these  solu¬ 
tions,  as  estimates  of  the  exact  solutions  but  with  the  reservations  mentioned  above.  It  should  be 
noted  that  in  both  cases,  the  original  authors  used  the  Richardson  extrapolation  method  as  out¬ 
lined  in  section  4  of  the  present  paper  to  estimate  some  of  the  uncertainties  in  their  results,  so  they 
may  be  considered  reliable. 
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2  NUMERICAL  METHOD 


The  equations  governing  the  steady,  two-dimensional,  incompressible  flow  and  heat  transfer  can 
be  written  in  dimensionless  variables  as: 


du  dv  . 

(1) 

du  du  dp  1 

(2) 

dv  dp  1  (d%  d\ )  T 

ay  dy  9^2  j  Fr 

(3) 

dT  dT  _  1  fa¥  aY) 

(4) 

where  Re,  Pr,  and  Fr  are  respectively,  the  Reynolds  number,  Prandtl  number  and  the  densimetric 
Froude  number. 

The  equations  are  solved  with  modified  versions  of  the  popular  TEACH  computer  code,  which  is 
based  on  the  SIMPLE  algorithm  of  Patankar  and  Spalding  (1972).  In  the  original  TEACH  code, 
the  differential  equations  were  discretized  using  numerical  approximations  to  derivatives  over  a 
staggered  computational  grid.  Convection  terms  were  approximated  with  a  hybrid  of  central  and 
first-order  upwind  differences,  depending  on  whether  the  cell  Peclet  number  was  less  or  greater 
than  two,  and  diffusion  and  other  terms  were  approximated  with  central  differences.  It  turns  out 
that  in  high  Reynolds  number  computations  on  typical  grids,  the  hybrid  scheme  mostly  degener¬ 
ates  to  a  first-order  upwind  scheme  which  introduces  “artificial  diffusion”  into  the  solution.  Effect 
of  the  discretization  of  the  convection  terms  is  considered  in  section  3,  where  the  hybrid  method  is 
compared  with  other  higher-order  difference  methods. 

The  boundary  conditions  for  case  1  and  2  are  shown  in  Figure  1  and  2.  For  the  inflow,  a  parabolic 
profile  is  prescribed  such  that  « (y)  =  \ly  -  24y^  for  x  =  0,0  <  y  <  112.  At  solid  surfaces,  the  typi¬ 
cal  “no  slip”  condition  of  zero  velocity  is  applied.  For  test  case  2,  additional  boundary  conditions 
are  necessary  for  the  temperature  equation.  The  non-dimensional  temperature  of  the  top  wall  is  set 
to  1 ,  while  that  of  the  bottom  wall  is  set  to  0.  A  linear  temperature  profile  is  prescribed  at  the  inlet 
plane  with  T  =  0  at  the  step  corner,  and  T  =  1  at  the  top  wall,  i.e.  Tiy)  =  2y  (x  =  0,  0  <  y  <  112). 
The  streamwise  temperature  derivative  is  set  to  zero  on  the  step  side  wall  (jc  =  0,  -ll2  <  y  <  0). 
The  boundary  condition  at  the  outflow  plane  will  be  described  in  section  5. 

In  sections  3  and  6,  uniform  grid  spacing  is  used  with  258  points  in  the  streamwise  direction  and 
34  points  in  the  cross  stream  direction.  An  x  domain  length  of  10  units  is  used.  For  the  results  of 
section  4  and  5,  nonuniform  grid  spacings  are  used.  Section  5  investigates  the  effect  of  the  loca¬ 
tion  of  the  outflow  boundary  by  using  four  x  domain  lengths.  Grid  1  is  defined  using  a  streamwise 
domain  length  of  30  units  and  480  streamwise  points.  Grids  2-4  are  simple  truncations  of  grid  1 . 
The  domain  lengths  for  grids  2-4  are  15,  10,  and  7  streamwise  units  giving  412,  372,  and  337 
streamwise  points,  respectively.  In  the  cross  channel  direction,  41  points  are  used  for  grids  1-4. 
The  streamwise  spacing  of  the  base  grid  1  is  expanded  geometrically  by  a  factor  of  1.01  going 
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from  0  to  30  units  downstream,  while  the  vertical  grid  spacing  is  expanded  by  a  factor  of  1 .05 
moving  from  the  top  and  bottom  walls  towards  the  center  line  of  the  channel.  Convergence  for  the 
model  problems  was  typically  achieved  in  about  4000  iterations. 

3  TRUNCATION  ERROR  IN  NUMERICAL  SCHEMES 

3.1  Description  of  the  numerical  schemes 

Truncation  error  in  a  numerical  scheme  may  result  from  errors  in  approximating  the  convection  or 
diffusion  terms,  but  in  high  Reynolds  number  flows  convection  usually  dominates  diffusion  so 
more  attention  needs  to  be  placed  on  the  former.  Further,  stability  or  algorithmic  considerations 
may  place  restrictions  on  the  form  of  the  convective  terms,  so  the  analysis  of  truncation  errors  in 
the  numerical  schemes  concentrates  on  the  approximation  for  the  convection  terms.  The  diffusion 
terms  are  simply  approximated  with  central  differences. 

Four  differencing  schemes  are  applied  to  the  two  model  problems.  The  differencing  schemes  are 
the  hybrid,  central,  second-order  upwind  (2nd  OU),  and  third-order  upwind  (3rd  OU).  On  dis¬ 
cretizing  the  governing  equations  (1),  (2),  (3)  and  (4)  over  a  typical  control  volume,  the  four 
schemes  lead  to  algebraic  equations  with  the  general  form: 

where  <|>  is  any  field  property  (velocity  or  temperature),  5u  is  the  source  term,  (pressure  gradient  or 
additional  viscous  terms),  the  A’s  are  the  convection/diffusion  coefficients,  and  B  is  the  fluid  body 
force,  where  applicable. 

A  typical  control  volume  is  shown  in  Figure  4  along  with  the  neighboring  points.  The  points  e,  n, 
w,  and  s  refer  to  values  at  the  east,  north,  west,  and  south  cell  faces,  respectively.  The  points  E,  W, 
N,  and  5  refer  to  grid  points  to  the  east,  north,  west,  and  south  of  the  cell  center,  respectively.  The 
points  EE,  WW,  MV,  and  SS  are  grid  points  located  to  the  east,  north,  west,  and  south,  two  points 
away  from  the  cell  center,  respectively. 

The  convection  through  the  cell  faces  is  defined  as: 


F,=  H,Ay 

(6) 

(7) 

Fn  =  v„A  X 

(8) 

(9) 

The  diffusion  at  the  cell  faces  is  defined  as: 

(10) 
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(11) 


3.1.1  Hybrid  scheme 

The  coefficients  for  the  hybrid  scheme  have  the  following  form: 


=  max  1^0,  D,  -  +  max  l-F^,  0] 

(14) 

=  max  1^0,  -  ^l^j]  +  max  [F^,  0] 

(15) 

=  max  [O,  D„  -  1|F„|]  +  max  [-F„,  0] 

(16) 

=  max  j^0,D^-  i|FjJ  +  max  [F,,  0] 

(17) 

^EE  ~  ^WW  -  ^NN  ~  ^SS  ~ 

(18) 

II 

a. 

(19) 

where  the  sum  for  Ap  is  taken  over  the  A  coefficients,  and  1 1  represents  the  absolute  value. 


A  truncation  of  the  hybrid  scheme  gives  the  first-order  upwind  (1st  OU)  scheme,  which  is  much 
less  accurate.  For  this  scheme,  the  coefficients  are  calculated  from  only  the  second  term  of  the 
right  hand  side  of  equations  (14)-(17).  If  the  cell  Peclet  number  was  greater  than  two  in  all  regions 
of  the  flow,  this  scheme  would  be  used  throughout,  but  such  a  condition  cannot  be  true  in  a  sepa¬ 
rated  flow. 
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3.1.2  Central,  second-order  upwind,  and  third-order  upwind  difference  schemes 
The  coefficients  for  the  higher-order  schemes  have  the  following  general  form: 


F, 

Ae  =  (D^  ~  y)  +  (/itnax  [F,,  0]  +/,max  [-F^,  0]  H-Zjmax  [-F,,  0] )  (20) 

=  (0.+  y )  +  Cfimax  [F,,  0]  +/,max  [-F^,  0]  -h/^max  [F^,  0] )  (21) 

F„ 

01  +/iinax  [-F^,  0]  -t-Z^max  [F„,  0] )  (22) 

F, 

=  (0,+  y)  +  </,max  [F„,  0]  +/imax  [-F„  0]  +f^max  [F,,  0] )  (23) 

^EE  =  -/ifnax  t"F,.0]  (24) 

^ww  ~  "/iniax  (F^,  0]  (25) 

^NN  =  -/i^ax  [-F„,0]  (26) 

= -/,max(F,.0]  (27) 

>4p  =  £>1,  (28) 


The  values  of/y  and/2  depend  on  the  higher-order  scheme  to  be  used  as  shown  in  Table  1. 


Table  l:Constants  used  in  A  coefficients  of  higher-order  schemes. 


Difference 

Scheme 

fl 

fl 

Central 

0 

0 

2ndOU 

0.5 

1.0 

3rdOU 

0.125 

0.25 

In  equations  (20)-(28),  the  second-  and  third-order  upwind  schemes  contain  two  groups  of  terms. 
The  first  group  contains  the  central  difference  terms,  while  the  second  group  contains  upwind- 
weighted  corrections  to  these  terms.  The  corrections  are  introduced  to  enhance  the  stability  of  the 
numerical  scheme.  Implementation  of  the  central  and  hybrid  schemes  requires  a  five-point  com¬ 
putational  stencil  at  each  node  whereas  the  others  require  a  nine-point  stencil.  To  minimize  algo¬ 
rithmic  changes,  all  the  higher-order  methods  are  implemented  via  deferred  correction  (see 
Khosla  and  Rubin  1974).  In  this  procedure,  the  coefficients  are  calculated  initially  using  equations 
(14)-(19)  for  the  hybrid  scheme.  As  the  solution  proceeds,  the  higher-order  scheme  is  slowly 
introduced  via  corrections  to  the  source  terms.  At  the  end  when  the  solution  is  fully  converged, 
the  coefficients  are  effectively  those  of  the  higher-order  scheme  outlined  in  equations  (20)-(28). 
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3.2  Results  for  different  approximations  of  the  convection  terms 


The  four  approximations  of  the  convection  terms  outlined  above  are  applied  to  the  two  test  prob¬ 
lems  over  a  streamwise  computational  domain  length  of  10.  Computed  locations  of  zero  shear 
stress  on  the  upper  and  lower  walls  are  shown  in  Table  2  for  case  1  and  Table  3  for  case  2.  From 
the  results  of  case  1,  the  hybrid  scheme  gives  a  large  percentage  difference  compared  to  the 
"benchmark"  result  (20%  for  the  top  point  and  17%  for  the  bottom  point).  The  artificial  diffusion 
introduced  by  the  hybrid  coefficients  makes  the  effective  Reynolds  number  lower  and  thus  the 
eddy  lengths  are  shorter  as  would  be  the  case  in  a  flow  with  slightly  lower  Reynolds  number.  To 
obtain  improved  accuracy  a  higher-order  method  must  be  used.  The  central  and  third-order 
upwind  schemes  give  percentage  differences  under  5%,  and  the  second-order  upwind  scheme 
under  8%. 

Similar  results  can  be  seen  for  the  points  of  zero  shear  stress  for  case  2.  The  hybrid  scheme  yields 
a  maximum  percentage  difference  (excluding  the  first  bottom  point)  of  around  25%.  The  second 
and  third-order  upwind  schemes  yield  percentage  differences  (excluding  the  first  bottom  point) 
under  5%  with  the  second-order  upwind  scheme  giving  the  closest  agreement  with  the  “bench¬ 
mark”  solution.  The  central  difference  scheme  yields  a  maximum  percentage  difference  under 
8%.  It  is  not  clear  why  the  second-order  upwind  results  agree  best  with  the  benchmark,  in  this 
case,  which  was  obtained  with  the  finite  element  method  with  bilinear  elements.  A  possible  expla¬ 
nation  would  be  that  both  numerical  schemes  are  roughly  similar  so  that  the  comparisons  contain 
a  measure  of  the  bias  in  the  numerical  schemes  rather  than  providing  a  pure  measure  of  accuracy. 


Table  2:  Points  of  zero  shear  stress,  case  1 


Difference 

Scheme 

Top(%  diff®) 

Bot.(%  diff®) 

Hybrid 

3.87(20) 

5.05(17) 

2ndOU 

4.47(7.8) 

5.70(6.6) 

Central 

4.64(4.3) 

5.88(3.6) 

3rdOU 

4.61(4.9) 

5.84(4.3) 

Benchmark 

4.85 

6.10 

a.  Percent  difference  between  value  and  benchmark  solution. 
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Table  3:  Points  of  zero  shear  stress,  case  2 


Difference 

Scheme 

lstTop(%'‘) 

2nd  Top 

3rd  Top 

1st  Bot*’ 

2nd  Bot 

3rd  Bot 

4th  Bot 

Hybrid 

1.01(17) 

4.38(25) 

-(-) 

0.65(45) 

2.14(17) 

3.76(21) 

6.70(25) 

2nd  OU 

1.17(3.3) 

5.50(.18) 

8.20(.49) 

0.36(0.0) 

2.43(2.9) 

4.54(.22) 

8.37(.24) 

Central 

1.14(6.1) 

5.23(5.0) 

7.65(6.7) 

0.44(18) 

2.41(3.7) 

4.32(5.3) 

7.99(4.9) 

3rd  OU 

1.16(4.3) 

5.32(3.2) 

7.78(4.9) 

0.43(16) 

2.43(2.9) 

4.39(3.6) 

8.10(3.6) 

Benchmaric 

1.21 

5.49 

8.16 

0.36 

2.50 

4.55 

8.39 

a.  Percent  difference  between  value  and  benchmark  solution. 

b.  This  comer  eddy  is  too  small  to  be  determined  accurately  on  the  scales  of  interest  here. 


The  order  of  the  difference  schemes  can  be  estimated  following  generalizations  of  the  Richardson 
extrapolation  method  to  be  described  in  detail  in  the  next  section.  Basically,  the  exact  functional 
value  can  be  approximated  in  terms  of  results  on  finite  grids  plus  the  leading  term  of  the  trunca¬ 
tion  error  as: 


<t»  =  + 

(29) 

(1.  =  02,+  (2/i)\  +  ... 

(30) 

<!'  =  <('4/.+  (4/1)  X+- 

(31) 

where  h  is  the  grid  spacing  in  the  jc-direction,  n  is  the  order  of  the  scheme  and  x„  is  a  grid  function, 
which  is  assumed  to  be  equal  for  the  h,  2h,  and  4h  grids.  The  grid  function  contains  spatial  deriv¬ 
atives  of  0  with  respect  to  x  and  y,  which  are  also  of  order  n.  The  statements  above  will  be  valid  so 
long  as  h  is  sufficiently  small  for  the  leading  term  to  be  dominant.  The  order  of  the  numerical 
scheme  can  then  be  estimated  from: 


n 


In 


In  (2) 


(32) 


Computations  were  made  on  the  standard  grid  (258  x  34)  and  two  sets  of  coarser  grids  (130  x  18 
and  66  x  10)  for  the  isothermal  backward  facing  step  (equivalent  to  case  1)  but  at  a  Reynolds 
number  of  100.  The  lower  Reynolds  number  was  chosen  for  economy  because  the  value  of  h 
required  for  equations  (29)-(32)  to  be  valid  is  larger  at  lower  Reynolds  numbers.  For  example,  the 
analysis  was  applied  to  computations  on  these  3  sets  of  grids  for  case  1  at  a  Reynolds  number  of 
400  using  the  hybrid  scheme.  Equation  (32)  estimated  the  hybrid  scheme  to  be  of  0th  order,  which 
is  clearly  incorrect.  Obviously,  at  this  Reynolds  number  the  coarser  grids  are  not  sufficiently  fine 
for  the  leading  term  in  the  truncation  series  to  be  dominant.  For  the  lower  Reynolds  number,  esti¬ 
mates  of  the  order  of  the  numerical  scheme  based  on  equation  (32)  are  presented  in  Table  4.  The 
results  show  that  the  1st  OU  scheme  is  indeed  first-order  accurate,  the  2nd  OU  and  central 
schemes  are  second-order  accurate,  and  the  3rd  OU  is  only  slightly  better  than  the  second-order 
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accurate.  The  only  possible  surprise  is  that  the  hybrid  scheme  is  almost  second-order  accurate. 
This  is  mainly  due  to  the  lower  Reynolds  number  of  the  flow  which  enables  the  use  of  central  dif¬ 
ferencing  in  major  sections  of  the  computational  domain.  The  ability  of  the  scheme  to  adapt  to 
flow  conditions  to  improve  its  order  of  accuracy  is  one  of  the  reasons  for  its  popularity,  but  the 
uncertainty  concerning  its  order  of  accuracy  in  different  flow  situations  makes  it  unreliable  for  use 
in  a  general  purpose  computer  code. 


Table  4:  Estimated  order  of  numerical  schemes,  case  \@Re=  100 


Difference 

Scheme 

1st  OU 

Hybrid 

2nd  OU 

Central 

3rdOU 

Order 

0.8 

1.9 

2.0 

2.0 

2.2 

4  DISCRETIZATION  ERROR 
4.1  Richardson  Extrapolation 

The  discretization  error  was  investigated  by  employing  a  Richardson  extrapolation  method  of 
analysis.  Grid  3  described  in  section  2  with  a  domain  length  of  10  units  was  used  as  a  base  grid. 
This  grid  will  be  defined  as  the  2h  grid.  A  grid  was  generated  with  half  the  grid  spacing  in  the  x 
and  y  directions  which  is  defined  as  the  h  grid,  and  one  was  generated  with  double  the  grid  spac¬ 
ing,  defined  as  the  4h  grid.  The  sizes  of  these  grids  are  displayed  in  Table  5. 


Table  5:  Grids  used  in  Richardson  extrapolation 


Grid 

Grid  size 

h 

743x81 

2h 

372x41 

4h 

187x21 

Central  differencing  was  used  for  the  convection  terms  of  the  governing  equations.  This  scheme, 
described  in  detail  in  the  previous  section,  is  second-order  accurate.  The  error  for  the  three  grids 
can  be  written  as  a  Taylor  series  following  Ferziger  (1989): 


=  h^X2  +  h^x^  +  ... 

(33) 

2  4 

^2h  ~  4h  X2+  16/l  X4+  ... 

(34) 

=  16/i^X2  + 256/i‘*X4+ ... 

(35) 
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If  h  is  sufficiently  small  so  that  the  first  term  on  the  right  hand  side  of  equations  (33)-(35)  domi¬ 
nates,  any  property  <t>  can  be  written  as: 


<t>  =  <!)*  +  £, 

(36) 

*1*  =  *t>2h  +  e2h 

(37) 

^  =  <t>4.  +  e4. 

(38) 

Equations  (33)-(35)  are  substituted  into  equations  (36)-(38)  respectively,  to  give: 


<t»  =  (t>^  +  h^X2 

(39) 

‘t’  =  <t'2.  + 

(40) 

^  =  ^4.  +  16/I^JC3 

(41) 

Equation  (39)  is  set  equal  to  (40): 

*t'.-<t>2.  = 

(42) 

From  equation  (33): 

e.  = 

(43) 

Equation  (43)  can  be  used  in  equation  (42)  to  obtain  an  expression  for  the  error  on  the  h  grid: 


The  error  estimate  given  by  equation  (44)  is  used  to  eliminate  the  h^X2  truncation  error  term  in 
equation  (33)  and  give  a  result  which  is  fourth-order  accurate  as: 


Note  that  this  result  is  generated  on  the  2h  grid,  since  equation  (45)  requires  a  fine  grid  {h)  and  a 
standard  grid  {2h)  value  at  each  location  the  equation  is  applied. 

Similarly,  one  can  combine  results  on  the  2h  and  4h  grids  to  yield: 


*1*4. 


4,. 


2nd 


2nd 


(46) 


Equation  (46)  gives  fourth-order  accurate  results  from  the  2h  and  4h  grid  values.  Sixth-order 
accurate  results  can  be  obtained  by  eliminating  the  second  term  in  equations  (33)-(35)  which 
assumes: 
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(47) 


‘'2* 


Aih 


Aih 


The  sixth-order  accurate  result  takes  the  form: 


hth 


nh\ 


*’■'''1  6lh 


I6h‘*x^ 

256/1^4 

I  + 

1 4/A  ^ 

+  256h^x, 

4/A  •* 


(48) 


(49) 

(50) 


By  equating  equations  (49)  and  (50)  and  following  the  procedure  outlined  above,  the  fourth-order 
accurate  truncation  error  is  estimated  as: 


<t'2A-‘!>4A 

15 


(51) 


The  fourth-order  accurate  result  can  be  corrected  to  obtain  a  sixth-order  accurate  result  as: 


f4A 


6th 


15‘*’2a|4,;, 


4th 


(52) 


Again  this  result  is  obtained  on  the  4h  grid.  The  procedure  used  to  generate  the  more  accurate 
results  in  equations  (45)  and  (52)  are  equivalent  to  those  for  estimating  zero-grid-size  results  (see 
Churchill  et  al.,  1981).  The  present  results  may  also  be  interpreted  as  the  Oh  results  for  the  second- 
order  scheme. 


4.2  Application  to  case  1 

As  an  example  of  the  estimation  of  discretization  error,  the  points  of  zero  shear  stress  at  the  upper 
and  lower  walls  are  calculated.  The  results  of  these  calculations  as  well  as  the  extrapolated  values 
using  equations  (45),  (46),  and  (52)  are  displayed  in  Table  6,  and  compared  to  the  "benchmark" 
solution. 
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Table  6:  Points  of  zero  shear  stress,  case  1 


Grid 

Top(%  diff^) 

Bot.(%  diff**) 

4.72(-2.6) 

6.01  (-1.5) 

4.55(-6.1) 

5.79(-5.1) 

3.61  (-25) 

4.76(-22) 

4.78(-I.4) 

6.08(-0.26) 

4^1  4,* 

4.87(0.41) 

6.13(0.52) 

4.77(-1.6) 

6.08(-0.31) 

Benchmark 

4.85 

6.10 

a.  Percent  difference  between  value  and  benchmark  solution. 

b.  Results  incompletely  converged  to  r,^  =  1.5xl0'^(maxi- 
mum  of  u,  V,  and  mass  equation  residual),  all  other  results 
converged  to  =  2.0x10  '*. 

The  results  using  Richardson  extrapolation  show  that  combining  the  lh\^^^ZinA  results, 

yields  an  improved  result.  For  the  top  wall  point,  the  percentage  difference  improves  from  6.1% 
(2/t  grid)  and  25%(4h  grid)  to  0.41%(4/i  grid,  fourth-order  accurate).  The  bottom  wall  point  per¬ 
centage  difference  improves  from  5.1%  {2h  grid)  and  22%{4h  grid)  to  026%(4h  grid,  fourth- 
order  accurate).  The  initial  results  on  the  h  grid  are  not  converged  to  the  same  level  as  the  2h  and 
4h  grid  results.  This  may  partly  explain  why  the  2h\  results  are  not  much  better  than  the  4/i|  4,,, 
results. 

4.3  Application  to  case2 

Similar  analysis  is  applied  to  test  case  2  with  calculations  completed  on  the  2h  and  4h  grid  only. 
The  h  grid  results  did  not  converge  fully  after  16,000  iterations  so  they  are  not  included  in  the 
analysis.  The  trends  are  similar  to  those  for  case  1.  The  results  are  shown  in  Table  7. 


Table  7:  Points  of  zero  shear  stress,  case  2 


Grid** 

lstTop(%^) 

2nd  Top 

3rd  Top 

IstBot*^ 

2nd  Bot 

3rd  Bot 

4th  Bot 

2^1  2nd 

1.19(1.7) 

5.30(3.5) 

7.75(5.0) 

0.47(31) 

2.46(1.6) 

4.39(3.5) 

8.10(3.5) 

4>«l  2nd 

1.02(16) 

4.75(16) 

6.93(18) 

0.62(42) 

2.21(13) 

3.92(16) 

7.45(11) 

1.25(3.2) 

5.48(.18) 

8.02(1.8) 

0.42(14) 

2.54(1.6) 

4.55(0.0) 

8.32(.84) 

Benchmark 

1.21 

5.49 

8.16 

0.36 

2.50 

4.55 

8.39 

a.  All  results  converged  to  =  5.0x10'^. 

b.  Percent  difference  between  value  and  benchmark  solution. 

c.  This  comer  eddy  is  too  small  to  be  determined  accurately  on  the  scales  of  interest  here. 
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5  OUTFLOW  BOUNDARY  CONDITIONS 


5.1  Description  of  boundary  conditions 

Most  of  the  boundary  conditions  are  straightforward  to  specify  as  discussed  in  section  2,  except 
for  the  location  of  the  outflow  plane.  The  usual  practice  is  to  locate  this  plane  far  enough  away 
from  the  region  of  interest  which  presumes  a  pre-knowledge  of  the  solution.  The  question  arises 
as  to  how  far  the  outflow  plane  should  be  located  in  separated  flows  and  what  errors  are  intro¬ 
duced  by  too  short  a  location.  The  effect  of  the  location  of  the  outflow  boundary  has  been  investi¬ 
gated  for  the  two  test  cases. 

The  outflow  boundary  condition  (OBC)  for  test  case  1  consists  of  setting  the  first  derivatives  of  u 
and  V,  in  the  direction  normal  to  the  outflow  plane,  to  zero,  while  satisfying  global  conservation  of 
mass  at  the  outflow.  Thus: 


^  "  j-0.5  ^"2.  j  -  1.;)  (53) 

^  =  (54) 


dv 

dx 


".7 


V 


n-  1.7 


(55) 


where  the  subscripts  2  and  n  denote  the  inflow  and  outflow  locations  respectively  for  the  u  vari¬ 
able.  If  global  continuity  is  satisfied  at  the  outflow  A  m  =  0  and  u„  j  =  y 

The  OBC  for  test  case  2  is  similar  to  that  of  test  case  1 ,  with  the  additional  boundary  condition  for 
the  temperature  equation  being: 


—  =  0  r  =  T 


(56) 


The  OBCs  for  the  two  model  problems  are  tested  by  obtaining  solutions  on  four  streamwise 
domain  lengths.  As  discussed  in  section  2,  the  four  domain  lengths  are  30, 15, 10,  and  7. 


5.2  Results  -  Case  I 


The  points  of  zero  wall  shear  stress  are  calculated  and  shown  in  Table  8,  while  the  streamlines  and 
pressure  contours  are  shown  in  Figures  5  and  6  for  the  various  domain  lengths.  It  is  clear  that  with 
the  outflow  boundary  conditions  specified  in  equations  (53)  to  (56)  the  location  of  the  outflow 
boundary  has  little  effect  on  the  computed  solution.  There  was  no  difficulty  in  obtaining  the  cor¬ 
rect  results  even  though  a  recirculating  eddy  was  dissected  by  this  boundary  for  the  domain  of  7. 
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Table  8:  Points  of  zero  wail  shear  stress,  case  1 


Domain 

Length 

30 

15 

10 

7 

%  diff^ 

1st  Top 

4.57 

4.56 

4.55 

4.53 

-0.66 

2nd  Top 

10.27 

10.27 

- 

- 

- 

1  St  Bot. 

5.80 

5.80 

5.78 

5.76 

-0.69 

a.  Percent  difference  between  long  (L=15)  and  short  (L=7)  domain. 


5.3  Results  -  Case  2 


The  points  of  zero  wall  shear  stress  for  case  2  are  shown  in  Table  9,  while  the  streamlines,  pres¬ 
sure,  and  temperature  contours  are  shown  in  Figures  7-9.  The  observation  above  of  little  influ¬ 
ence  of  the  location  of  outflow  boundary,  is  also  true.  As  shown  by  Wfllson  et  al.  (1991),  other 
boundary  conditions  such  as  zero  second  derivatives  normal  to  the  outflow  plane  are  not  as  suc¬ 
cessful,  Only  the  present  boundary  conditions  (equations  53-56)  enabled  the  pressure  variation 
across  the  outflow  plane  in  Figures  6  and  8  to  be  correctly  predicted  for  the  shortest  domain  length 
of  7. 


Table  9:  Points  of  zero  wall  shear  stress,  case2 


Domain 

Length 

30 

15 

10 

7 

%diff® 

1  St  Top 

1.18 

1.18 

1.19 

1.20 

1.7 

2nd  Top 

5.29 

5.29 

5.30 

5.24 

-0.95 

3rd  Top 

7.91 

7.91 

7.75 

- 

- 

4th  Top 

10.19 

10.19 

- 

- 

- 

1st  Bot. 

0.47 

0.47 

0.47 

0.47 

0.0 

2nd  Bot 

2.45 

2.45 

2.46 

2.47 

0.82 

3rd  Bot. 

4.38 

4.38 

4.39 

4.35 

-0.68 

4th  Bot. 

8.25 

8.25 

8.10 

- 

- 

a.  Percent  difference  between  long  (L=I5)  and  short  (L=7)  domain. 


13 


6  INCOMPLETE  ITERATIVE  CONVERGENCE 


6.1  Stopping  criteria  and  error  estimation 

Uncertainty  due  to  incomplete  iterative  convergence  can  be  defined  as  the  difference  between  the 
current  and  the  exact  solution  of  the  discretized  problem  on  the  same  grid.  The  discretized  solu¬ 
tion  will  never  satisfy  the  continuous  equations  exactly.  Therefore,  there  will  be  a  point  in  the  iter¬ 
ation  process  when  further  relaxation  of  the  system  of  equations  will  not  bring  any  additional 
improvement  in  the  solution.  Stopping  criteria  must  be  selected,  and  the  resulting  errors  from 
foregoing  additional  iterations,  must  be  estimated.  One  stopping  criterion,  which  is  common  in 
practice,  is  to  terminate  the  iterative  process  when  the  difference  in  computed  results,  from  one 
iteration  to  the  next,  falls  below  a  pre-selected  amount.  Another  stopping  criterion  is  based  on  a 
measure  of  how  well  the  discretized  solution  satisfies  the  discretized  equations.  This  quantity  is 
referred  to  as  the  residual  of  the  discretized  equation.  In  this  method,  when  the  residual  falls 
below  a  pre-selected  amount,  the  iterations  are  stopped. 

A  method  has  been  proposed  by  Ferziger  (1989)  that  is  based  on  not  only  the  difference  of  succes¬ 
sive  iterates  but  also  on  the  rate  of  convergence.  The  uncertainty  due  to  incomplete  iterative  con¬ 
vergence  can  also  be  estimated.  The  solution  after  the  iteration  can  be  written  as: 


(t)"  =  4»  +  e'‘  (57) 

where  <t)"  is  the  solution  after  n  iterations,  4)  is  the  desired  converged  solution,  and  e"  is  the  error 
due  to  incomplete  convergence.  The  iterative  scheme  for  a  linear  system  can  be  simplified  as: 

(58) 

where  A  is  the  amplification  matrix  that  transforms  the  solution  after  n  iterations,  4)",  to  (t)"*  *  after 
adding  an  effective  source  term,  S.  If  equations  (57)  and  (58)  are  combined,  it  can  be  shown  that 
the  error  obeys  the  homogeneous  form  of  equation  (58).  The  eigenvalues  and  eigenfunctions  are 
calculated  from: 


The  initial  error  e°  can  be  written  as  a  generalized  Fourier  series: 


*  =  I 

Equation  (60)  can  be  used  in  (57)  to  develop  an  expression  for  the  solution  after  n  iterations  as: 

n 

4.'- =  0+ (61) 

*=  1 

If  it  is  assumed  that  X,  will  dominate  after  many  iterations,  equation  (61)  can  be  replaced  with  a 
simplified  expression: 
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<t)"  =  <}i  +  a,A.;y,  (62) 

Equation  (62)  can  be  written  using  the  indices  n-1,  n,  and  n+1  to  yield  an  estimate  for  the  princi¬ 
pal  eigenvalue; 


a"*  ‘  A 

<l>  -4> 

a"  a"  “  1 

<t>  -<)> 


An  estimate  of  the  principal  eigenvalue,  applicable  to  all  nodes,  can  be  obtained  as  the  ratio  of  the 
L2-norms  of  the  numerator  and  denominator  of  equation  (63).  If  equation  (62)  is  rearranged,  an 
estimate  of  the  convergence  error  can  be  calculated  as: 


The  result  of  equation  (64)  is  that  the  convergence  error  depends  on  the  difference  in  the  solution 
from  iteration  to  iteration  plus  the  rate  of  convergence.  This  means  that  a  given  difference 
between  successive  iterates  will  yield  different  convergence  errors  depending  on  the  principal 
eigenvalue  for  the  iteration  matrix.  Therefore,  if  the  iterative  method  exhibits  poor  convergence, 
the  principle  eigenvalue  of  the  iteration  matrix  A  will  be  close  to  unity.  This  will  make  it  difficult 
to  reduce  the  convergence  error  e"  because  the  inverse  of  the  denominator  of  equation  (64)  will  be 
large.  This  result  can  form  the  basis  for  stopping  criteria.  When  the  convergence  error,  defined  by 
equation  (64),  falls  below  a  predetermined  level,  the  iterations  can  be  terminated. 


6.2  Application  to  model  problems 

The  convergence  error  was  estimated  for  cases  1  and  2  by  using  the  u  component  of  velocity  for 
the  general  variable  <ti  in  equations  (57)-(64).  The  convergence  error  is  shown  in  Figure  10.  In 
addition,  the  norm  of  the  residuals  of  the  u  momentum,  v  momentum,  and  continuity  equations  is 
shown  in  Figure  1 1 .  The  evolutions  of  the  recirculation  zones,  determined  by  the  points  of  zero 
shear  stress  at  the  upper  and  lower  walls,  are  also  monitored  throughout  the  iterative  process. 
They  are  shown  in  Figure  12  and  13  for  cases  1  and  2,  respectively.  This  gives  an  indication  when 
additional  computational  work  will  not  change  the  solution. 

In  Figure  12,  the  first  point  of  zero  wall  shear  stress  for  case  1  on  the  top  wall  and  the  length  of  the 
bottom  recirculation  zone  remain  constant  after  approximately  3000  iterations.  This  implies  that 
the  shear  stress  along  the  top  wall  is  relatively  constant  and  that  the  solution  is  not  changing.  From 
Figure  11,  the  norm  of  the  u,  v,  and  mass  equation  residuals  at  this  point  is  2  x  10'-^  and  the  esti¬ 
mated  convergence  error  is  about  1  x  10"^.  It  appears  that  this  is  an  appropriate  stopping  point. 

As  shown  in  Figure  13,  the  three  recirculation  zone  lengths  for  case  2  remain  constant  after 
approximately  3200  iterations.  Note  that  the  second  bottom  eddy  is  the  last  to  attain  a  constant 
length.  Figure  1 1  shows  that  the  norm  of  the  residuals  at  this  point  is  approximately  6  x  10'^,  and 
the  estimated  convergence  error  is  approximately  6  x  10'\ 

Figures  10  and  1 1  have  similar  trends  and  values.  For  case  1,  the  norm  of  the  residuals  (Figure  1 1 ) 
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follows  roughly  the  same  constant  slope  throughout  the  iterative  process.  The  mean  value  of  the 
convergence  error,  (Figure  10)  also  follows  the  same  slope  shown  in  Figure  11.  The  values  of  the 
norm  and  convergence  error  are  also  similar.  For  case  2,  the  residual  norm  follows  the  same  slope 
as  case  1  until  about  800  iterations,  after  which  it  decreases  at  a  slower  rate.  The  mean  value  of 
the  convergence  error  also  shows  the  same  behavior.  The  residuals  shown  in  Figure  1 1  have  been 
plotted  using  the  L2-norm  of  the  u  momentum,  v  momentum,  and  continuity  equations,  while  Fig¬ 
ure  10  uses  the  error  estimate  of  the  u  component  of  velocity. 

For  the  two  test  cases,  both  the  norm  of  the  residuals  and  the  estimate  of  the  convergence  error 
appear  to  be  appropriate  stopping  criteria. 

7  COMPUTATIONAL  GRID  ASPECT  RATIO 


7.1  Choice  of  grid-cell  aspect  ratio 

In  most  computational  fluid  flow  problems,  the  choice  of  the  grid-cell  aspect  ratio  is  not  trivial. 
Convergence  characteristics  suggest  that  the  aspect  ratio  should  be  of  order  unity,  but  the  need  to 
resolve  boundary  layers  may  dictate  much  higher  aspect  ratios.  The  natural  way  to  carry  out  grid 
refinement  is  to  halve  the  cell  size  in  each  direction,  which  automatically  maintains  the  initial  cell 
aspect  ratio.  What  is  the  effect  of  this  choice  on  the  accuracy  of  computed  results? 

Assuming  that  the  leading  truncation  error  term  is  of  second-order  and  that  the  grid  spacing  is  suf¬ 
ficiently  small,  the  functional  value,  <|),  can  be  written  as: 


hld\  h]d\ 

where  and  hy  are  the  grid  spacing  in  the  x  and  y  directions,  respectively. 

Let  the  aspect  ratio  be  defined  as,  /!/;  =  hjhy  also  let  p  =  — /  —  then: 

h-f  dx 


The  discretization  error  is  then  approximately: 


For  the  same  total  number  of  grid  points,  refinement  may  be  selectively  in  the  x  or  y-direction. 
The  choice  would  produce  a  change  in  the  cell-aspect  ratio.  For  example,  if  hy  is  reduced  by  a  fac¬ 
tor  m,  with  /ijj  {=h)  unchanged,  Af^  will  be  increased  by  the  factor  m.  Then  the  discretization  error 
would  be  approximately: 
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On  the  other  hand,  if  hy  is  unchanged  while  is  reduced  by  a  factor  m  so  that  =  him  and  Aj^  is 
also  reduced  by  a  factor  m,  the  error  is: 


(69) 


The  ratio  of  the  errors  is  then; 


E. 


The  implication  of  this  result  is  that: 

^-KifV=P  (71) 

^<l,ifA/<  p 

‘"2 

^>i,ifV>  p 

Therefore  the  limiting  cell  aspect  ratio  for  selective  grid  refinement  is  or  So  long  as 

A/;  <  Vp  grid  refinement  in  the  y-direction  leads  to  more  effective  error  reduction  than  in  the  x- 
direction.  If  >  Vp  further  grid  refinement  in  the  y-direction  becomes  less  effective  than  that  in 
the  x-direction.  In  boundary-layer  type  flows  p  »  1,  so  the  optimum  value  of  Ag>  1.  But  in  sep¬ 
arated  flows  with  no  preferred  direction  p  ~  1  and  the  optimum  A/^  will  also  be  about  1 . 


1  + 


rnAji 
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(70) 


7.2  Application  to  model  problems 


To  explore  these  effects,  computations  of  case  1  were  made  on  several  grids  with  cell  aspect  ratios 
in  the  range  of  0.625  to  5.0.  Both  hybrid  and  central  difference  schemes  were  utilized  because  it 
was  difficult  to  get  converged  solutions  with  the  latter  on  coarser  grids  because  of  wiggles  gener¬ 
ated  by  the  well-known  “odd-even”  decoupling  problem.  The  results  for  the  hybrid  scheme  are 
presented  in  Table  10  and  those  for  the  central  scheme  in  Table  11.  Grids  1-3  have  the  same  cell 
size  in  the  x-direction,  but  the  y-direction  cell-sizes  were  halved  in  2  and  halved  again  in  3.  Grids 
4  and  5  have  the  same  y-direction  cell-sizes  as  grids  2  and  3,  respectively,  but  doubled  x-direction 
cell-sizes.  Correspondingly,  the  aspect  ratios  are  doubled.  The  results  show  improved  agreement 
with  the  benchmark  solution  with  refinement  in  the  y-direction  corresponding  to  increased  aspect 
ratio.  For  example,  in  Table  10  grids  2  and  5  have  the  same  total  number  of  points  but  the  case 
with  the  higher  aspect  ratio  gives  better  agreement  with  the  benchmark.  This  confirms  the  analysis 
for  large  p.  The  same  applies  to  grids  1  and  4.  Results  in  Table  1 1  for  the  central  difference 
scheme  show  that  the  deviation  from  the  benchmark  was  reduced  by  a  factor  of  4  simply  by  halv¬ 
ing  the  cell-size  in  the  y-direction.  For  the  second-order  scheme,  such  a  reduction  would  normally 
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be  expected  from  halving  of  the  cell  size  in  both  directions.  Clearly,  the  y-component  of  the  trun¬ 
cation  error  is  dominant  and  selective  refinement  in  this  direction  is  more  cost  effective  than  a  glo¬ 
bal  refinement.  This  cannot  be  presented  as  a  panacea  for  all  two-dimensional  separated  flows,  but 
will  depend  on  the  characteristics  of  the  flow  and  the  extent  of  deviation  from  an  elongated 
boundary-layer  character.  The  closer  large  regions  of  the  flow  are  to  boundary-layer  type  flows, 
and  hence,  the  larger  p  is,  the  more  will  be  the  tendency  for  large  cell  aspect  ratios  to  produce 
more  accurate  solutions.  For  separated  flows  with  little  or  no  elongated  regions  with  boundary- 
layer  type  flows,  aspect  ratios  of  order  unity  would  be  the  most  effective. 


Table  10:  Points  of  zero  wall  shear  stress,  hybrid  difference  scheme,  case  1. 


Grid 

Grid  Size 

Aspect  Ratio 

1st  Top(%^) 

1st  Bot 

1 

258  X  18 

0.625 

2.36(51.3) 

3.52(42.3) 

2 

258  X  34 

1.25 

3.87(20.2) 

5.04(17.4) 

3 

258  X  66 

2.5 

4.53(6.6) 

5.76(5.6) 

4 

130x34 

2.5 

3.33(31.3) 

4.42(27.6) 

5 

130x66 

5 

4.12(15.1) 

5.30(13.1) 

Benchmark 

- 

- 

4.85 

6.10 

a.  Percent  difference  between  value  and  benchmark  solution 


Table  11:  Points  of  zero  wall  shear  stress,  central  difference  scheme,  case  1. 


Grid 

Grid  Size 

Aspect  Ratio 

lstTop(%^) 

1  St  Bot 

1 

258  X  18 

0.625 

„b 

__b 

2 

258  X  34 

1.25 

4.64(4.3) 

5.88(3.6) 

3 

258  X  66 

2.5 

4.80(1.0) 

6.05(0.8) 

Benchmark 

- 

- 

4.85 

6.10 

a.  Percent  difference  between  value  and  benchmark  solution 

b.  Solution  not  converged. 


8  CONCLUDING  REMARKS 

Various  sources  of  uncertainty  in  numerical  computations  of  fluid  flow  have  been  examined.  Spe¬ 
cific  estimates  of  numerical  error  magnitudes  were  computed  with  reference  to  two  two-dimen¬ 
sional  separated  flow  problems.  Truncation  error  in  numerical  schemes  can  be  estimated  by 
comparing  solutions  from  low  and  higher-order  schemes.  The  effect  of  outflow  boundary  condi¬ 
tions  can  be  estimated  by  varying  systematically  the  location  of  the  outflow  boundary  without 
changing  the  grid  distribution  or  the  numerical  scheme.  Discretization  errors  can  be  estimated  by 


making  computations  on  related  grids  with  varying  degrees  of  fineness  and  using  Richardson 
extrapolation  method.  The  solution  can  then  be  improved.  The  method  can  also  be  used  to  deter¬ 
mine  the  global  order  of  accuracy  of  a  numerical  method.  The  uncertainty  in  computed  results  due 
to  incomplete  convergence  of  the  iterative  scheme  can  be  removed  by  computing  an  estimate  of 
the  convergence  error  and  using  this  as  a  stopping  criterion  rather  than  the  more  widely  used 
change  in  computed  results  between  iterates.  Grid  aspect  ratio  effects  on  the  solution  are  also 
important.  Higher  aspect  ratios  are  more  effective  in  generating  accurate  solutions  in  separated 
flows  with  elongated  regions  with  boundary  layer  character. 
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Figure  1.  Geometry  for  case  1,  Isothermal  backward  facing  step. 


Figure  2.  Geometry  for  case  2,  Stratified  backward  facing  step.  Same  velocity 

boundary  conditions  as  case  1. 
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Figure  6.  Normalized  pressure  contours,  case  1 
a)  "Benchmark"  solution,  b)  L=7,  c)  L=10,  d)  L=15,  e)  L=30 

Level  values  are:  0.01,  0.02,  0.03,  0.04, 0.05, 0.06, 0.07, 0.08,  0.09,  0.10,  0.12, 0.14,  0.16,  0.01 
0.20, 0.22, 0.24. 
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Figure  8.  Normalized  pressure  contours,  case  2 
a)  "Benchmark"  solution,  b)  L=7,  c)  L=10,  d)  L=15,  e)  L=30 

Level  values  are:  0.00  at  step  corner,  0.04841  increments. 
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Iterations 


Figure  11.  L2-Norm  of  residuals  of  the  U  mom.,  V  mom. 
and  continuity  equations  for  model  problems. 
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